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Molecular dynamics simulations are carried out to investigate the diffusion behavior of penetrable- 
sphere model fluids characterized by a finite energy barrier e. The self-diffusion coefficient is eval- 
uated from the time-dependent velocity autocorrelation function and mean-square displacement. 
Detailed insights into the cluster formation for penetrable spheres are gained from the Enskog fac- 
tor, the effective particle volume fraction, the mean free path, and the collision frequency for both 
the soft-type penetrable and the hard-type reflective collisions. The simulation data are compared 
to theoretical predictions from the Boltzmann kinetic equation and from a simple extension to finite 
e of the Enskog prediction for impenetrable hard spheres (e — > oo). A reasonable agreement between 
theoretical and simulation results is found in the cases of e* = e/ksT = 0.2, 0.5, and 1.0. However, 
for dense systems (packing fraction <j> > 0.6) with a highly repulsive energy barrier (e* = 3.0), a 
poorer agreement was observed due to metastable static effects of clustering formation and dynamic 
effects of correlated collision processes among these cluster-forming particles. 

PACS numbers: 51.20.+d, 05.20.Dd, 61.20.Ja, 05.60.-k 
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I. INTRODUCTION 

The so-called penetrable-sphere (PS) pair potential is 
defined as 



(1) 



where a is the diameter of the penetrable spheres and e > 
is the strength of the repulsive energy barrier between 
two overlapping spheres when they penetrate each other. 
The PS model was suggested by Marquest and Witten [l[ 
as a simple theoretical approach to micelles in a solvent 
to explain the experimentally observed crystallization of 
copolymer mesophases, where a simple cubic solid phase 
coexisted with the disordered suspension. This simple 
model system has been the subject of several theoretical 
and simulation studies An excellent review in this 

area up to 2001 can be found in Ref. [ll|. More recently, 
the PS model has been extended to include a short-range 
attractive tail [THH]. 

The PS interaction model reduces to the classical hard- 
sphere (HS) system when e* = e/ksT — > oo (where 
T is the temperature and ks is the Boltzmann con- 
stant). This is equivalent to the zero-temperature limit 
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T* = ksT/e —> 0. In the opposite high-penetrability 
or high-temperature limit (e* — > 0, T* — > oo) the PS 
system becomes a collisionless ideal gas. Except in the 
pure HS case, penetrability allows one in principle to con- 
sider systems with any value of the nominal packing frac- 
tion <f> = (ir/6)na 3 , where n is the number density. The 
first correction to the equilibrium ideal-gas structural and 
thermodynamic properties in the combined limit e* —> 0, 
qy — > oo with e*qy = const has been exactly obtained Q. 

Recently, the nonequilibrium transport properties of 
unbounded potentials, such as the linear-core (l5j and 
Gaussian-core [l6| fluids have received much attention. 
In the case of the PS model, the Liouville operator and 
the Boltzmann-Lorentz collision operator were long time 
ago derived in Rcfs. 



171 ] and [18(, respectively. As a 
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more explicit application of kinetic theory to soft mat- 
ter one of the authors has derived the relevant 
transport coefficients (self-diffusion, shear viscosity, and 
thermal conductivity) in the context of the Chapman- 
Enskog method for the Boltzmann equation of dilute 
gases ((f) <C 1) [20j , Interestingly, in the PS binary col- 
lision dynamics the particle penetration is analogous to 
the double refraction of light through a sphere made of a 
transparent material of relative refraction index depend- 
ing on the relative collision velocity and the repulsive 
energy parameter [lj| . 

In addition to transport properties, two of us [2l| have 
applied two different theoretical predictions, based on 
the fundamental-measure theory proposed by Schmidt 
Q and the bridge density-functional approximation pro- 
posed by Zhou and Ruckenstein [22|, to the inhomoge- 
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neous structure of PS model fluids in the spherical pore 
system. More recently [23[ , as a continuation of theoreti- 
cal approaches along this direction, the modified density- 
functional theory (based on both the bridge density func- 
tional and the contact-value theorem) has been investi- 
gated for the structural properties of PS fluids near a slit 
hard wall and the Verlet-modified bridge function for one- 
component systems proposed by Choudhury and Ghosh 
Q has also been extended to PS fluid mixtures. 

Theoretical results for such a precisely defined model 
potential can be directly compared against the exact 
machine-experimental data obtained from computer sim- 
ulations via Monte Carlo (MC) and molecular dynamics 
(MD) methods. So far, almost all simulations for the PS 
model fluid have been carried out using the MC method. 
On the other hand, better statistics can be achieved in 
MD simulations, particularly in systems with discontin- 
uous interactions. For instance, in order to calculate the 
virial route to the equation of state for hard-core systems, 
MC computations require an accurate estimation of the 
radial distribution function at the contact point. Com- 
putationally, the pair distribution function may change 
rapidly near the contact distance in the systems of ionic 
solutions, highly charged colloids, aligned liquid crystals, 
etc. Under those conditions the extrapolation of the pair 
distribution function to the contact value may lead to 
large uncertainties. For this reason, the pressure deter- 
mined by MC simulations is known to be less accurate 
than that by MD computations [2~ij . 

The main motivation in the present work is to un- 
dertake a detailed molecular-based simulation study of 
transport properties of PS model fluids, which can be 
in turn directly compared with theoretical and/or em- 
pirical predictions. MD data for the PS interaction po- 
tential have not been presented in the literature before. 
More specifically, we have investigated the PS model sys- 
tem via the equilibrium MD method over a wide range 
of densities and repulsive energy parameters to investi- 
gate the equation of state, the self-diffusion coefficient, 
and its related time-dependent quantities including the 
velocity autocorrelation function (VACF) and the mean- 
square displacement (MSD). Together with existing ap- 
proximations for the equilibrium equation of state as well 
as for the self-diffusion coefficient in the Boltzmann and 
Enskog-likc descriptions, our simulation results can be 
used to assess and construct a fundamental basis of theo- 
retical and practical predictions for the relevant transport 
properties. Such simulation approaches at the atomic or 
molecular level can also be used in improved statistical 
integral-equation theories of liquid state and help in in- 
terpreting real experimental data. 

The organization of this paper is as follows. Section 
Hi! presents some simple theoretical approximations for 
the equation of state and the self-diffusion coefficient of 
the PS fluid. The MD method employed in this paper is 
briefly described in Sec. IIII1 The most important part of 
the paper is contained in Sec. IIV1 where the simulation 
results for the pressure, the self-diffusion coefficient, the 



Enskog x factor, the effective particle volume fraction, 
the mean free path, the collision frequencies, the veloc- 
ity autocorrelation function, and the mean-square dis- 
placement are presented, compared with theoretical ap- 
proaches, and discussed. The paper is closed with some 
concluding remarks in Sec. [Vj 



II. THEORETICAL APPROACHES FOR THE 
PS FLUID 

A. Equation of state 

The equilibrium compressibility factor Z = p/nk^T in 
the PS model, where p is the pressure, is given by 



^PS 



1 + 4(f>xx 



PS 



where 



(2) 



(3) 



is a parameter measuring the degree of penetrability of 
the particles (ranging from x = in the free-penetrability 
limit e* — > to x = 1 in the opposite impenetrability 
limit e* — > oo) and x = 5 PS (c + ) is the contact value 
of the radial distribution function g (r). 

In Ref. [lfl ] two approximate theories were proposed 
to obtain <? PS (r): one valid in the high-penetrability 
(i.e., small-e*) regime and the other one in the low- 
penetrability (i.e., large-e*) regime. We will refer to 
these two theories as the high-penetrability approxi- 
mation (HPA) and the low-penetrability approximation 
(LPA) , respectively. We give below the expressions for x 
in both approximations. 

In the HPA, x is given by [l(| 



x PS = i 



xw(l)e xw{1 \ 



(4) 



where 



480a; f°° „ (k cos k - sin k) sin(fcr) 
nr J k 3 — 240a; (k cos k — sin k) k 2 

(5) 

Equation Q reduces to the exact result x = 1 + xw(l) 
in the limit x — > (i.e., e* — > 0) with <f>x = const. 
In the case of the LPA, one has [l(| 
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Here, A is obtained from the transcendental equation 
(l-x)A 2 0(1 + 0/2) 
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FIG. 1: (Color online) Different possible trajectories in a 
collision with a scaled relative velocity u* 2 = 1.1. The 
dashed trajectories correspond to "soft" (refraction) colli- 
sions, while the solid trajectories correspond to "hard" (spec- 
ular) collisions. In general, the collisions are of hard type if 
v *2y/l ~ (b/cr) 2 < 1. As the (reduced) repulsive energy bar- 
rier e* increases, less and less collisions are soft. 
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1,2,3) are the three roots of the cubic 
Bi, B 2 , and B% are 
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In the limit x — > 1 (i.e., e* — > 00), the solution of Eq. 
([7| is A = 1, and so one recovers the solution of the 
Percus-Yevick integral equation for HS [25[ . 



B. Self-diffusion coefficient 

In the low-density regime (0 — > 0) the transport coeffi- 
cients of a gas made of particles interacting via a poten- 
tial u(r) can be derived by application of the Chapman- 
Enskog method to the well-known Boltzmann equation 
[2fjj . The crucial quantity distinguishing an interaction 
potential from another one is the scattering angle as a 
function of both the impact parameter b and the relative 
velocity of the colliding pair, v\ 2 . 

In the particular case of the PS potential [lj|, if the 
scaled relative velocity v* 2 = vi 2 / ye/m, where m is the 
mass of a particle, is smaller than 1, the "projectile" par- 
ticle does not have enough kinetic energy to penetrate 
into the core of the "target" particle, and consequently 
it is deflected exactly as if the target were a hard sphere. 

On the other hand, if v* 2 > 1 and b/a < \J 1 — v* 2 -2 , i.e., 

if u* 2 yl — (b/a) 2 > 1, the projectile traverses the core 



analogously to the double refraction of light through a 
sphere made of a transparent material of relative refrac- 
tion index \ 1 — v$ 



In fact, if b/a > J 1 — v* 



i.e., 1 < v* 2 < 1/yl — (b/a) 2 , a "total reflection" pro- 
cess takes place and the deflection is again as in the HS 
case. Figure [T] shows different possible trajectories corre- 
sponding to v* 2 = 1.1. 

In the first Soninc approximation, the self-diffusion co- 
efficient -Do obtained from the Boltzmann equation for 
the PS model is given by [l9j 
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where 
with 

Ri(y) 



m (j) fl^ 
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(y-l)(y + 2) , 4y 2 -4y + 3 



(9) 
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Obviously, in the low-penetrability limit (e* — > 00), 
the self-diffusion coefficient for the PS model in Eq. © 
reduces to that of the HS model, namely 



£n HS = ^ 



1 InkBTa 
16 ' 



(12) 



As said above, Eqs. © and (fT2"j) are derived from the 
Boltzmann equation (in the first Sonine approximation) , 
and thus they are well justified in the high-dilution limit 
0—^0. On the other hand, they do not account for finite- 
density effects. To correct this deficiency, several empir- 
ical or semi-empirical expressions have been proposed in 
the case of the HS model. Among them, the most basic 
one is provided by the Enskog kinetic theory [2(|. The 
Enskog correction for the self-diffusion coefficient in the 
HS system is 
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(13) 



where the Enskog factor x is the contact value of the 
radial distribution function of the HS fluid. This quantity 
is related to the corresponding equation of state in terms 
of the compressibility factor by 



Z HS = 1 



40X 



HS 



(14) 



Equation (|13|) takes into account that the effective num- 
ber of collisions in a dense gas is increased by a factor 
X HS . Consequently, the self-diffusion coefficient is de- 
creased by the same factor, relative to the Boltzmann 
prediction at the same density. An excellent approxi- 
mation for x HS within the stable fluid region (0 < < 
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0.494) is provided by the Carnahan-Starling (CS) for- 
mula [25L ual 



HS 



0/2 



(1 



(15) 



There are also a number of empirical formulas for the 
D HS . For systems of 500 HS particles or slightly fewer, 
the following analytical fit to MD data was reported by 
Speedy f27|: 



D 



HS 



D HS 



1 + Ci 



C2 



(16) 

Here, (j) g ~ 0.57 is the packing fraction at the HS glass 
transition and Speedy's values are c\ = 0.48 and c 2 = 
1.17. Recently, a much more extensive MD computation 
was performed by Sigurgeirsson and Heyes [28| with an 
efficient MD algorithm dealing with up to 32000 HS par- 
ticles. They refined the values of the fitting coefficients 
ci and c 2 in Eq. (pi]) as ci = 0.4740 and c 2 = 1.1657. 
The empirical form (|16[) takes into account the crowding 
effects in the first bracket term and the hydrodynamic 
back-flow effects at intermediate densities in the second 
bracket term. Both the Enskog formula (|13|) and the 
empirical expression (|16[) have in common the fact that, 
as expected on physical grounds, D HS < Z?^ s , i.e., the 
self-diffusion coefficient decays more rapidly than hypcr- 
bolically with increasing density. 

In the case of the PS system, the task of extending 
the Boltzmann result © to finite density to get the self- 
diffusion coefficient D PS is much more difficult than in 
the HS case. In fact, the ratio D ps /Dq S is not only a 
function of density (as in the HS case) but also a function 
of temperature or, cquivalently, of e*. Based on the En- 
skog result (TiUj) , it seems natural to propose the following 
Enskog-like expression: 



D 



PS 



nPS 

,,ps : 



(17) 



where Dq s is given by Eqs. (|^|)— (jTT|) and \ PS 1S cither 
obtained from the empirical values of Z PS from Eq. @ 
or derived from the HPA [Eqs. §4§ an d (J5J] r from the 
LPA [Eqs. ©-©]. According to Eq. (HTJ) , L» PS < D PS . 
However, as will be seen in Sec. IIV1 this inequality is in 
general not supported by our MD data. 



III. COMPUTATIONAL METHOD 

As a successful diagnostics tool, molecular-based com- 
puter simulations are usually employed to investigate the 
underlying diffusion behavior of the model system of in- 
terest. To this end, in this work we have carried out 
microcanonical MD simulations for the PS model fluid 
in a manner similar to that originally proposed by Alder 
and Wainwright for hard-core systems |29j , which is well 



described elsewhere [24| . Post-collision velocities for col- 
liding pairs of particles are assigned according to the type 
of collision (see Fig.Q}: either hard-type specular reflec- 
tion or soft-type refraction. In all cases both the total 
momentum and kinetic energy are conserved in these PS 
collision conditions. 

The initial configurations with 864 penetrable spheres 
were generated by randomly inserting particles with ve- 
locities drawn from Maxwell-Boltzmann distributions. 
The initial configurations were aged, or equilibrated, for 
5 x 10 7 collisions before accumulating the final simula- 
tion results. Additional ensemble averages were evalu- 
ated from a total number of 5 x 10 8 collisions. 

Our MD algorithm has been tested in a number of 
ways. When the repulsive energy parameter was assigned 
a large value (typically e* > 3) at the low-density regime 
{<p < 0.2), the static and dynamic results generated from 
our MD simulations faithfully reproduced the pure HS 
system. Our resulting MD calculations for a few selected 
runs were also compared with MC computations reported 
in the literature. A good agreement with previous MC 
data for the thermodynamic and structural properties 
[Io| again confirmed the validity of the MD method em- 
ployed in this work. All MD results reported here are 
scaled to dimcnsionlcss quantities by using a unit parti- 
cle diameter er, a unit particle mass m, and a unit ther- 
mal energy ksT. In these system units the reduced self- 
diffusion coefficient is expressed as D* = Dja^k^Tjm. 



IV. RESULTS AND DISCUSSION 
A. Equation of state 

Before presenting the results for the self-diffusion co- 
efficient, which is the main quantity of interest in this 
work, it is worth considering the equation of state. 

Figure [2] compares the MD simulation data for Z with 
the HPA and the LPA. For the most penetrable case 
(e* = 0.2) both theories agree well with the MD data, 
with the agreement being excellent in the case of the HPA 
(relative deviations between MD and HPA smaller than 
0.1%). For e* = 0.5 the HPA is still quite good, while for 
€* = 1 the LPA performs very well. It is remarkable that 
both theories, while being based on opposite approaches, 
are so close each other up to e* m 1 and densities as 
large as <j> = 1. In the least penetrable case (e* = 3) 
the LPA behaves reasonably well up to (f> = 0.3 (where 
the PS system is only slightly distinguishable from a HS 
system, represented here by the CS equation of state) 
but strongly overestimates the MD values for larger den- 
sities, when penetrability effects start to play a relevant 
role. By accident the HPA does a good job for e* = 3 
and <f> w 1. It is important to remark that the MD data 
for the PS fluid with e* = 3 clearly deviate from the HS 
values for cf> > 0.2. The reason is that, even though the 
value e* — 3 represents a rather high energy barrier, as 
the density increases more particles are forced to overlap, 
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FIG. 2: (Color online) Compressibility factor Z versus the 
packing fraction <j> for several values of e* . The circles are MD 
results, the solid lines are the HPA predictions [Eqs. ([2jl-(f5jl] 
and the dashed lines are the LPA predictions [Eqs. ((2)1 and 
©-(H}]. The diamonds at e* = 1 and & = 0.4, 0.5, 0.6, and 
0.8 are MC simulation data from Ref. For comparison, 

the curve corresponding to the HS fluid described by the CS 
equation of state [Eqs. (|14p and (|15p ] is also plotted (dash- 
dotted line). 



which results in a substantial decrease in the pressure rel- 
ative to that of the HS system at the same density. 



B. Self-diffusion coefficient 

The self-diffusion coefficient is a single-particle quan- 
tity which has been more frequently studied in MD cal- 
culations than other collective transport properties, such 
as the shear viscosity and the thermal conductivity. The 
self-diffusion coefficient D can be determined from the 
temporal integration of the VACF using the Green-Kubo 
formula 



D 



1 

3./, 



<a<v(o)-v(t)> 



(18) 



and also from the slope of the MSD versus time using the 
Einstein relation 



(r 2 (t)) =6Dt. 



(19) 



In our MD simulations these two methods have produced 
consistent results, typically with less than 3% differences. 

By using a semi-logarithmic scale, in Fig. [3] we display 
the reduced self-diffusion coefficient D* as a function of 
the packing fraction 0, as obtained from MD simulations, 
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FIG. 3: (Color online) Reduced self-diffusion coefficient D* 
as a function of the packing fraction <f> for several values of e* . 
The symbols are MD simulations for the PS system, the solid 
lines are the Boltzmann theoretical predictions in the high- 
dilution limit for the PS fluid [cf. Eqs. ©HUH]; the dotted 
lines are the Enskog-like theoretical approximations for the PS 
fluid [cf. Eq. I|17|l]. complemented with the empirical values of 
X PS ; the dash-dotted line (interrupted at the maximum pack- 
ing fraction of HS systems, <f> = 7rv2/6 — 0.7405) represents 
the Boltzmann theoretical values in the high-dilution limit for 
the HS fluid [cf. Eq. (|12|l ]; the dash-dot-dotted line represents 
the Enskog theoretical values for the HS fluid [cf. Eq. (|13|) ]. 
complemented with the CS values of x HS j an d t- ne thick solid 
line represents the empirical fitting-values from MD simula- 
tions for the HS system [cf. Eq. (Q2J with ci = 0.4740 and 
c 2 = 1.1657]. 



from the Boltzmann theoretical predictions in the high- 
dilution limit [Eqs. (!§))- (ITl"]) ] and from the Enskog-like 
theoretical approximations [Eq. (|17p]. For comparison, 
Fig. [3] also includes curves corresponding to the HS sys- 
tem (e* — > oo ): the Boltzmann theoretical values in the 
high-dilution limit [Eq. (|12p], the Enskog theoretical val- 
ues [Eq. (|T3| ]. and the empirical fitting- values from MD 
simulations [Eq. (flB]) ] with c\ = 0.4740 and c 2 = 1.1657. 
There are several interesting features that can be ob- 
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served in Fig. [3] As expected, the self-diffusion coeffi- 
cient decreases with increasing density. It also decreases 
with increasing repulsive energy barrier at fixed <f>. This 
behavior is not counterintuitive. As the (reduced) en- 
ergy barrier e* increases, more and more collisions are 
of hard (specular) type. This gives rise to an increase 
in the effective collision frequency and thus a decrease 
in diffusion. At the Boltzmann level of description, this 
effect can be grasped from comparison between Eqs. (j9]) 
and (fT2]l and the property Jlfj < f . For given values of 
e* and <f>, the self-diffusion coefficient evaluated from the 
Boltzmann kinetic equation is larger than the one from 
the Enskog-type correction equation. As shown by Eq. 
(fTF]) . both coefficients are simply related by the \ factor, 
which is the contact value of the radial distribution func- 
tion, i.e., x = g( (J+ )- This value is close to 1.0 in the 
high-dilution limit and increases with increasing density 
due to particle crowding effects near the contact distance. 

Except for the case e* = 3.0, the qualitative behavior of 
the MD self-diffusion data is similar for the different val- 
ues of e*. In the cases e* = 0.2, 0.5, and 1.0, we observe 
in Fig.[3]that the PS Boltzmann approximation produces 
a reasonably good agreement with the corresponding MD 
data, even for relatively large values of the packing frac- 
tion (where the Boltzmann approximation would be ex- 
pected to break down). In the case e* = 0.2 the relative 
deviations of the Boltzmann predictions with respect to 
the MD data increase with density. For e* = 0.5, how- 
ever, the relative deviations reach a maximum at about 
<p = 0.6 and then slowly decay with increasing density. 
In the two previous cases the Boltzmann prediction [Eqs. 
(|9|)- (|11|) ] underestimates the self-diffusion coefficient, at 
least in the density domain < <f> < 1. This qualitative 
feature changes in the case e* = 1.0, where the Boltz- 
mann values are below the MD ones up to <f> — 0.5 only. 
In fact, the best general agreement between the Boltz- 
mann results and the simulation data occurs, because of 
this accidental crossing, for e* = 1.0. It is expected that, 
on increasing the barrier, the crossing shifts toward lower 
densities. For the highly repulsive system with e* = 3.0, 
there is no crossing effect, and so the Boltzmann approxi- 
mation overestimates the MD values of D* , being reliable 
only in the narrow range of densities (f> < 0.2. 

For the pure HS fluid, it has been reported [H!,[3(| that 
reliable self-diffusion data, significantly improving over 
the Boltzmann values, are obtained from the Enskog ki- 
netic equation within the range of equilibrium stable HS 
fluids (0 < 0.494). One may see this point from Fig. [3] 
by comparing the MD-fitting diffusion data with the HS 
Boltzmann and Enskog theoretical approximations. It 
would then seem natural to expect a similar improve- 
ment of the Enskog-like correction (JTTJ) over the bare 
Boltzmann prediction also in the case of the PS fluid. In- 
terestingly enough, however, this expectation only turns 
true for the highest repulsive barrier considered (e* = 3, 
where a good agreement with MD data is observed in 
the density range (f) < 0.5), as well as for e* = 1 and 
4> > 0.8. Since, according to Eq. (JTTJ) , the Enskog-likc 



values of D* arc smaller than those obtained from the 
Boltzmann equation, the former values are worse than 
the Boltzmann values whenever the latter underestimate 
the MD data. As shown in Fig. El this is what happens 
for e* = 0.2 and 0.5, as well as for e* = 1 and <f> 0.5. 

The system with e* = 3.0 deserves further comments. 
First, as indicated above, neither the Boltzmann equa- 
tion nor the Enskog-likc correction provides reliable val- 
ues of D* for cj) > 0.6. In fact, the MD values decay with 
increasing density much more rapidly than both theories 
predict. Moreover, the MD diffusion data (in the semi- 
logarithmic representation in Fig. E|) are seen to exhibit 
an inflection point near (f> = 0.4, a qualitative feature not 
accounted for by either the Boltzmann or the Enskog the- 
ories of the PS fluid. On the other hand, the existence of 
an inflection point of log D* vs <f> is also present in the HS 
system, as observed from the Enskog and the empirical 
lines in Fig. El In fact, the MD data of the PS fluid with 
e* = 3 are close to the HS values up to <\> s=s 0.2, analo- 
gous to what happens with the pressure (see Fig. [5]). For 
larger densities, however, the self-diffusion of PS parti- 
cles is much larger than that of HS particles at the same 
packing fraction. One may suggest that, in this range 
of higher densities with a high repulsive interaction, the 
self-diffusion process is greatly affected by static struc- 
tural effects together with dynamic correlated motion in- 
volved in the PS model system. Later in this section, 
we will comment on some interesting observations from 
additional structural and dynamic results obtained from 
our MD calculations. It may be worthwhile noting here 
that, even in the simple HS system, there are various 
transition properties, mostly investigated by simulation 
approaches, such as freezing of the fluid ((f) = 0.494), 
melting of the solid ((f> = 0.545), loose random packing 
((f) = 0.555), glass transition (<f) = 0.57), and dense ran- 
dom packing (0 = 0.64) [H. 



C. Enskog x factor and thermodynamic consistency 

In order to gain more detailed insights into the thermo- 
dynamic and structural properties related to the diffusion 
behavior, we have also examined other relevant proper- 
ties of the PS model interaction system. Here, we start 
by considering the Enskog-type correction x factor. 

According to statistical mechanics, there are two main 
routes in theoretical and simulation studies to obtain \ in 
a PS system: (a) one from the contact value of the radial 
distribution function as x PS = g PS (a + ) and (b) another 
one from the pressure via the so-called virial route as rep- 
resented by Eq. i.e., x PS = (Z PS - 1)/40(1 - e" e *). 
In simulation approaches, x values in the contact-value 
method (a) can be directly computed during simulation 
runs, while in the second method (b) they are indirectly 
evaluated from MC ensemble averages or MD time av- 
erages of the compressibility factor at a given thermo- 
dynamic condition. In principle, except for unavoidable 
computational errors, those two methods should generate 
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FIG. 4: (Color online) Enskog-correction factor xasa func- 
tion of packing fraction <f> f° r several values of e* , as obtained 
in MD simulations. The solid symbols are obtained from the 
contact-value method, while the open symbols are obtained 
from the compressibility factor. The solid line represents the 
function ~x HB {4 1 ) given by the CS equation of state for the HS 
fluid [Eq. J15}]. 



the same value in equilibrium stable liquid states. 

As shown in Fig. 21 an excellent agreement between 
those two methods is observed over the entire density 
range for e* = 0.2, 0.5, and 1.0. Again, the system 
with e* = 3.0 requires separate comments. In that case, 
we observe that, up to ~ 0.2 the MD values are very 
close to the HS contact values obtained from the CS for- 
mula. This agrees with what is observed in Figs. [5] and [3] 
Nevertheless, the HS values of x increase rapidly as the 
density increases, while the PS values reach a maximum 
around <f> = 0.4 and decrease thereafter. The most strik- 
ing feature observed in Fig. @] is the separation between 
the MD values of x obtained from methods (a) and (b) 
with e* = 3.0 and <fi > 0.6, with the maximum relative 
deviation being of almost 30% at (j> = 1. 

Also, although not shown in Fig. @J we have evaluated 
the values of the ratio g(a + ) /[g(a~)e' i ] by using the ex- 
trapolated MD data from g(r). This ratio should take the 
unity value, except for computational errors, at any given 
density and repulsive energy parameter if the system is 
in an equilibrium stable liquid state. We have observed 
that the deviations of g(a + )/[g(a~)e e } from unity in the 
cases e* = 0.2 and e* = 0.5 are less than 0.1% and 0.3%, 
respectively. The internal agreement between g{<J + ) and 
g(<r~)e e is also very good in the case e* = 1.0, with 



FIG. 5: (Color online) (a) Effective particle volume fraction 
4> e s and (b) ratio 4> c s/<f> as functions of the packing fraction 
(j>. The symbols (with lines to guide the eye) represent MD 
data for several values of e*. The upper and lower dashed 
lines correspond to the limiting cases of the non-overlapping 
HS system (e* — > oo) and the totally random overlapping PS 
system (e* — ¥ 0), respectively. 



a maximum deviation of about 3% at (f> = 0.9. Again, 
the least penetrable case (e* = 3.0) presents a peculiar 
behavior. Up to </> = 0.5 the ratio g(<r + ) /\g{a~)e t ) devi- 
ates from 1 less than 5%, but thereafter it markedly in- 
creases with density until having g(a + )/[g(a~)e e ] ~ 1.6 
at <f> = 1. 

All these discrepancies between the values of x ob- 
tained from g{a + ) 1 (Z - l)/4</>(l - e~ e *), and g{o~)e e * 
when e* = 3.0 and </> > 0.6 are likely due to the ap- 
pearance of cluster- forming metastable structures. Under 
these conditions we have employed the direct contact- 
value method x = <7( <7+ ) to obtain the Enskog-type 
kinetic approximations for the self-diffusion coefficient 
shown in Fig. [3] 



D. Effective particle volume fraction 

As a dimensionlcss measure of the number of particles 
per unit volume we are using the nominal packing frac- 
tion cj> = (ir/6)ncr 3 . In the case of HS systems, c/> coincides 
with the fraction of the total volume that is occupied by 
the spheres. For PS systems, however, particles can inter- 
penetrate, thus reducing the fraction of occupied volume. 
Let us define the effective particle volume fraction <fi e ff as 
the average effective total volume occupied by PS par- 
ticles divided by the system volume. Of course, (freff is 
a function of both <fi and e* that satisfies the inequality 
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<j>cft < with the equality taking place only in the HS 
limit (e* — > oo) or in the low-density limit (<f> — > 0). For 
fully penetrable systems in the high-penetrability limit 
(e* — > 0), the corresponding particle configuration will 
become that of a totally unbiased random structure, and 
this leads statistically to cf) e g = 1 — 

We have calculated </> e ff by the simple hit-and-miss 
method using a uniform (10<7 x lOcr x lOcr) grid over 
approximately half a million equilibrium configurations 
during our MD computations. The results are displayed 
in Fig. [5l For the systems with e* = 0.2, 0.5, and 1.0, 
the resulting data are very close to (but slightly above) 
the curves representing randomly distributed configura- 
tions. In the case e* = 3.0, the MD values of 4> e g are 
close to 4> up to 4> « 0.2, which indicates HS-like config- 
urations in that density range. As the density increases 
further, 4> e s/4> significantly decreases and at <fi « 0.6 it 
even crosses the random distribution expectation. It is 
paradoxical that at <f> = 1, for instance, particles are so 
much overlapped in the case e* = 3.0 that less than 55% 
of the available volume is actually occupied by them; in 
contrast, at the same density, particles occupy about 65% 
of the total volume in the case of a much less repulsive 
barrier e* = 0.2. Again, this signals in the case e* = 3.0 
a transition near <j> = 0.6 to a metastable state charac- 
terized by a high degree of clustering. 



E. Mean free path and collision frequencies 
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FIG. 6: (Color online) (a) Reduced mean free path A* — X/a 
(in semi-logarithmic scale) and (b) product \*<j> as functions 
of the packing fraction. The symbols (with lines to guide the 
eye) represent MD data for several values of e*. The dashed 
and dash-dotted lines correspond to the HS system in the 
Boltzmann and Enskog approximations, respectively. 



In the HS kinetic theory, the mean free path A HS in 
the dense system is related to that in the high-dilution 
limit A^ s = 0-/6^0 by the x HS factor in a similar way 
as in Eq. @3), i.e., A HS = A[p/x HS [H- In this case 
of HS systems the mean free path simply characterizes 
the typical distance (much higher than the diameter a 
in the dilute regime) traversed by a sphere between two 
successive collisions. Of course, the distance between the 
centers of the two particles in a HS collision is r = a + . In 
contrast, two classes of events contribute to the mean free 
path in the PS system (see Fig. [J). On one hand, one still 
has hard collisions taking place at r = <r + . On the other 
hand, soft encounters give rise to a primary "external" 
collision at r = a + followed by a secondary "internal" 
collision at r = cr~; this second event contributes to the 
mean free path with a value smaller than a. Therefore, 
every soft collision has two contributions (primary and 
secondary) to the mean free path. For dilute systems, 
where A > a, the primary contribution is on the order 
of A, while the secondary one is practically zero; as pen- 
etrability increases, almost all the collisions are soft and 
thus the mean free path is about half the value obtained 
in a HS system at the same (low) density. 

We have evaluated the mean free path for PS fluids 
during MD simulations, and the results are displayed in 
Fig. |U In the semi- logarithmic scale of Fig. [HJa), the 
values of A* = X/a display a similar decaying behavior 
for the cases e* = 0.2, 0.5, and 1.0. However, again in 



the case e* = 3.0 a different behavior can be observed, 
where A* is seen to exhibit a weak density dependence for 
4> > 0.6. Figurc[6jb) shows that, in the case e* = 3.0, the 
product X*(j> agrees with the HS prediction up to <j> = 0.2, 
but then it presents an accentuated minimum at (f> ~ 
0.4; this peculiar behavior is another distinct evidence 
for cluster-forming structures. On the other hand, the 
curves for e* = 0.2, 0.5, and 1.0 clearly depart from the 
HS values, even for small densities. This illustrates the 
penetrability effects in the PS fluid and the influence of 
soft collisions, as discussed above. For instance, in the 
case e* = 0.2 the values of X*4> slowly decay with density, 
taking a nearly constant value that is almost half the 
Boltzmann HS value (\^ s /a)4> = 1/6^ ~ 0.11785. 

For HS systems, the collision frequency is w HS = 
2(v)/X HS , where (v) = ^8fc s T/ 7rm is the average speed. 
In the case of PS systems, it is instructive to decompose 
the collision frequency u> = uj s + lj^ into the soft-type 
collision frequency (u> s ) and the hard- type collision fre- 
quency (u>h)- They are calculated as follows. Every time 
a collision with relative speed v±2 and impact parameter b 
occurs, it is cataloged as either soft or hard depending on 
whether (mt)j 2 /e)[l — (b/a) 2 ] — 1 is positive or negative, 
respectively. If J\f s (t) and Afh (t) denote the total num- 
bers of soft and hard collisions, respectively, over a time 
interval t, then w s = 2J\f s (t)/tN and ui h = 2J\f h (t)/tN, 
where N is the total number of particles. 

The collision frequencies uj s and uih are plotted in Figs. 
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FIG. 7: (Color online) Reduced soft-penetrable collision fre- 
quency lu* = u) s g I y/ksT '/m and (b) reduced hard-reflective 
collision frequency oj^ = u)hO / ' \JksT '/in as functions of the 
packing fraction <f>. The symbols (with lines to guide the eye) 
represent MD data for several values of e*. The dashed line 
in (b) corresponds to the HS system, i.e., oj* hs = (A8/y/n)(j). 



Eta) andEJb), respectively. We observe that soft colli- 
sions are dominant in the cases e* = 0.2, 0.5, and (to 
a lesser extent) 1.0. In particular, soft collisions repre- 
sent about 90% of all collisions for all the systems with 
e* = 0.2. It is interesting to note that both uj s and toh 
grow almost linearly with <fr in those three cases. In con- 
trast, if e* = 3.0 we find an initial quasi-linear growth 
followed by a much slower regime for > 0.6. As might 
be expected, for e* = 3.0 hard collisions dominate over 
soft ones, with the former representing about 90% of all 
collisions. 



F. Velocity autocorrelation function and mean 
square displacement 

Both the VACF [see Eq. flT5Jl] and the MSD [sec 
Eq. (fT9]) ] provide useful insights into the dynamic time- 
dependent behavior related to diffusion processes in PS 
systems. It is then illustrative to analyze the manner 
in which those functions change with increasing densities 
and repulsive energy parameters. 

In Figs. [Sfa) and MX>) we display the VACF (in units 
of 3fcgT/m) as a function of time in units of the relevant 
HS mean collision time t hs = 1/cj hs = A HS /2(u) for 
e* = 0.5 and e* = 3.0, respectively, and several charac- 
teristic densities. The primary mechanism involved in the 
rapid decay of the VACF is provided by hard collisions, 
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FIG. 8: (Color online) Normalized VACF as a function of the 
reduced time in units of r HS for (a) e* = 0.5 and (b) e* = 3.0. 
The curves correspond to cf> = 0.1 (solid lines), 0.3 (dashed 
lines), 0.5 (dash-dotted lines), 0.7 (dash-dot-dot lines), and 
0.9 (dotted lines). Note the different horizontal scales in both 
panels. 



so that colliding particles rapidly forget their initial ve- 
locities through successive collisions. For soft-penetrable 
collisions, post-collision velocities (or, equivalently, the 
colliding particle trajectories) are relatively correlated 
with their own initial values. In agreement with this, 
the normalized VACFs for the systems with e* = 0.5 ex- 
hibit similar exponentially decaying behaviors for differ- 
ent densities, since the soft-penetrable collision process 
is dominant in this case. As a consequence the areas be- 
low the curves are hardly dependent on <f>, which implies 
D ~ t hs oc (jr 1 . In the case e* = 3.0 the resulting VACF 
for (f> = 0.1 exhibits a decaying exponential behavior sim- 
ilar to that of the system with e* = 0.5 and the same den- 
sity, except that now the decay is much more rapid due 
to the prevalence of hard collisions. For higher densities, 
the development of cluster- forming structure significantly 
influences the collective motion of penetrable particles by 
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FIG. 9: (Color online) Log-log plot of the MSD (in units of 
a 2 ) as a function of the reduced time in units of t hs for (a) 
e* = 0.5 and (b) e* = 3.0. The curves correspond to <p — 0.1 
(solid lines), 0.3 (dashed lines), 0.5 (dash-dotted lines), 0.7 
(dash-dot-dot lines), and 0.9 (dotted lines). 



retardation (for incoming particles) or acceleration (for 
outgoing particles) of the velocities of colliding pairs. Un- 
der those conditions the resulting VACFs exhibit a non- 
exponential behavior. Furthermore, the trajectories of 
colliding particles are largely restricted by backscatter- 
ing or cage effects between clusters in the higher-density 
regime {4> > 0.6 and e* = 3.0), as shown in Fig.(5Jb). 

As displayed in the log-log plot of Fig. HI trends sim- 
ilar to the ones mentioned above for the VACF can be 
observed from the MD results for the MSD curves. At 
very short times, before hardly occurring particle col- 
lision, the MSD curves changes more rapidly than in 
longer times. This is due to the ballistic motion of par- 
ticles until they collide with their neighbors, resulting in 
(r 2 (t)) ~ f 2 . Subsequent collisions make trajectories re- 
semble a random-walk diffusion process, so that, after a 
certain transition period, the diffusive regime defined by 



Eq. (fT9|) is established. These behaviors are clearly illus- 
trated in Figs.^a) and[^b), except that the MSD curve 
for <j) = 0.9 and e* = 3.0 presents a transient anomalous 
sub- diffusive behavior where (r 2 (t)) ~ £ 7 , with 7 < 1. 
This signals a hindrance of the diffusion process by ob- 
struction or trapping phenomena. Once the diffusive 
regime (7 = 1) is reached for longer times, it is character- 
ized by a very low value of the self-diffusion coefficient, 
as shown in Fig. [3] 



G. Discussion on the Boltzmann kinetic theory 

Before concluding this section, it is of interest to return 
to one of the observations made earlier to explain some 
relevant shortcomings involved in the Boltzmann theo- 
retical approximation. As illustrated in Fig. [31 together 
with Figs. H] and |H1 the failure of the Boltzmann kinetic 
approximation for the PS model fluid becomes more im- 
portant as the density increases. This is not surprising: 
one may recall that the Boltzmann kinetic theory is based 
on the high-dilution limit. A key element related to this 
kinetic theory is the molecular chaos assumption, known 
as "Stosszahlansatz," in which the pre-collision velocities 
of two colliding particles are assumed to be totally uncor- 
rected. In addition, regardless of a given model poten- 
tial, the Boltzmann kinetic theory deals with only binary 
collision effects by totally neglecting multiple collisions. 
As observed in MD simulations for the PS model poten- 
tial in this work, the deviation between our MD diffusion 
data and the Boltzmann predictions can be largely due 
to the neglect of such spatio-temporal correlations in the 
PS collision dynamics, particularly in dense system with 
cluster-forming structures. For example, a simple conjec- 
tural argument will intuitively give, in analogy with the 
HS relation £> HS cx A , that particles with larger mean 
free paths become more diffusively dispersed (larger diffu- 
sion coefficients), and vice versa. However, for the dense 
PS fluid with the highest repulsive barrier {<p > 0.6 and 
e* = 3.0), our MD results clearly manifest contradic- 
tions against this simple conjecture: similar values of the 
mean free path (Fig. [5]) yield a significant reduction in 
the corresponding self-diffusion coefficient (see Fig. [3]). 
This partly explains that the Boltzmann kinetic theory 
does not take proper account of the consequences of cor- 
related collisions on the self-diffusion coefficient in the PS 
system, as investigated in this work. More detailed MD 
simulation studies can be very helpful to enable quali- 
tative predictions of the underlying behavior of the PS 
interaction systems together with statistical-mechanical 
approaches, which will be one of our current research 
topics for further theoretical and simulation work. 



V. CONCLUDING REMARKS 

In this work, as an intermediate between theory and 
experiment, MD simulations have been carried out to 
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investigate the detailed diffusion behavior of PS model 
fluids. The self-diffusion coefficient has been calculated 
from its related time-dependent properties of the VACF 
and the MSD. The resulting simulation data have been 
used to assess theoretical predictions by the Boltzmann 
kinetic equation and an Enskog-likc correction. Detailed 
insights involved in the cluster formation for penetrable 
spheres have been observed from the effective particle 
volume fraction, the mean free path, and the collision 
frequency for both the soft-type penetrable and the hard- 
type reflective collisions. 

A reasonable good agreement with Boltzmann and En- 
skog theoretical approximations is found in the cases 
e* = 0.2, 0.5, and 1.0. On the other hand, for the 
dense PS fluid with the highest repulsive barrier (</> > 0.6 
and e* = 3.0), several distinct evidences of the cluster- 
forming structure are exhibited from static structural 
and dynamic collisional properties. In that case, a poor 
agreement between theory and simulation is observed due 
to those effects, especially correlated collision processes. 
Under those conditions, the indirect virial route and the 
direct contact- value route to obtain the Enskog-type cor- 



rection factors yield inconsistent results. This indicates 
that for e* = 3.0 and 4> > 0.6 the states reached in our 
MD simulations are not of strict thermodynamic equilib- 
rium but are long-lived metastable states. We are cur- 
rently examining such cluster-formation conditions with 
larger numbers of penetrable spheres to check the num- 
ber dependence of transport properties including the self- 
diffusion, the shear viscosity, and the thermal conductiv- 
ity coefficients. 



Acknowledgments 

S.-H.S. wishes to express his thanks to Jae-Moon Yang, 
Young- Jin Ha, and Jong-Hoon Sohn for their help and 
implementations during simulation runs. The research 
of A.S. was supported by the Ministerio de Educacion y 
Ciencia (Spain) through Grant No. FIS2007-60977 (par- 
tially financed by FEDER funds) and by the Junta dc 
Extremadura through Grant No. GRU10158. 



C. Marquest and T. A. Witten, J. Phys. France 50, 1267 [19 
(1989). 

C. N. Likos, M. Watzlawek, and H. Lowen, Phys. Rev. E 
58, 3135 (1998). 

M. Schmidt, J. Phys.: Condens. Matter 11, 10163 (1999). [20 
M. J. Fernaud, E. Lomba, and L. L. Lee, J. Chem. Phys. 
112, 810 (2000). 

M. Schmidt and M. Fuchs, J. Chem. Phys. 117, 6308 [21 
(2002). 

N. Choudhury and S. K. Ghosh, J. Chem. Phys. 119, [22 
4827 (2003). 

L. Acedo and A. Santos, Phys. Lett. A 323, 427 (2004). [23 
Al. Malijevsky and A. Santos, J. Chem. Phys. 124, 
074508 (2006). [24 
A. Santos and Al. Malijevsky, Phys. Rev. E 75, 021201 
(2007); 75, 049901(E) (2007). [25 
Al. Malijevsky, S. B. Yuste and A. Santos, Phys. Rev. E 

76, 021504 (2007). [26 
C. N. Likos, Phys. Rep. 348, 267 (2001). 
A. Santos, R. Fantoni, and A. Giacometti, Phys. Rev. E [27 

77, 051206 (2008). [28 
R. Fantoni, A. Giacometti, Al. Malijevsky, and A. San- 
tos, J. Chem. Phys. 131, 124106 (2009). [29 
R. Fantoni, A. Giacometti, Al. Malijevsky, and A. San- 
tos, J. Chem. Phys. 133, 024101 (2010). [30 
A. Donev, B. J. Alder, and A. L. Garcia, J. Stat. Mech.: 
Theory Exp. (2009) P11008. [31 
L. A. Shall and S. A. Egorov, J. Chem. Phys. 132, 184504 
(2010). [32 
A. R. Altenberger, Physica A 80, 46 (1975). 

L. Groome, J. W. Dufty, and M. J. Lindenfeld, Phys. 
Rev. A 19, 304 (1979). 



A. Santos, in Rarefied Gas Dynamics: 24th International 
Symposium on Rarefied Gas Dynamics, AIP Conf. Proc. 
No. 762, edited by M. Capitelli (AIP, Melville, NY, 2005), 
pp. 276-281. 

S. Chapman, T.G. Cowling, The Mathematical Theory of 
Nonuniform Gases (Cambridge University Press, Cam- 
bridge, England, 1970). 

S.-C. Kim and S.-H. Suh, J. Chem. Phys. 117, 9880 

(2002) . 

S. Zhou and E. Ruckenstein, J. Chem. Phys. 112, 8079 
(2000). 

S.-C. Kim, B.-S. Seong, and S.-H. Suh, J. Chem. Phys. 
131, 134701 (2009). 

M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Clarendon, Oxford, 1987). 

J.-P. Hansen and I. R. McDonald Theory of Simple Liq- 
uids (Academic Press, Amsterdam, 2006). 
N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 
635 (1969). 

R. J. Speedy, Mol. Phys. 62, 509 (1987). 

H. Sigurgeirsson and D. M. Heyes, Mol. Phys. 101, 469 

(2003) . 

B. J. Alder and T. E. Wainwright, J. Chem. Phys. 31, 
459 (1959). 

B. J. Alder, D. M. Gass, and T. E. Wainwright, J. Chem. 
Phys. 53, 3813 (1970). 

S.-H. Suh, W.-K. Min, and J. M. D. MacElroy, Bull. 
Korean Chem. Soc. 20, 1521 (1999). 
P. Cutchis, H. van Beijeren, J. R. Dorfman, and E. A. 
Mason, Am. J. Phys. 45, 970 (1977). 



